RNAvelocity4：velocyto.R的使用

分享是一种态度

RNA注释

这里以inDrop实验数据举例，spliced/unspliced的RNA可以通过：

1. 使用dropEst 输出管道生成 类似10X的 bam 文件：

~/dropEst/build/dropest -m -F -L eiEIBA -o run1 -g cellranger/refdata-cellranger-mm10-1.2.0/genes/genes.gtf -c ~/dropEst/configs/indrop_v3.xml *.bam

1. 使用velocyto.py来注释spliced/unspliced的reads，写出标准loom文件：

velocyto run -u Gene -o out -e SCG71 -m mm10_rmsk_srt.gtf -v SCG_71_tophat.filtered.sorted.bam UCSC/mm10/Annotation/Genes/genes.gtf

（注意，也可以使用-V参数直接通过 DropEst 注释spliced/unspliced的reads：）

~/dropEst/dropest -V -C 6000 -m -g ucsc_mm10_exons.gtf.gz -c ~/dropEst/configs/indrop_v3.xml *.aligned.bam)

下面的例子从 velocyto.py 生成的loom文件开始，使用 pagoda2^[1] 获取细胞clusters/embedding，然后估计/可视化RNA速率。

加载R包：

library(velocyto.R)

Loading required package: Matrix

数据下载

下载预先计算的loom矩阵

wget http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom

读取矩阵

ldat <-read.loom.matrices(url("http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom"))
#Error: is.character(name) is not TRUE

使用pagoda2标准化和聚类细胞

使用spliced 表达矩阵作为pagoda2的输入，并查看分布。

emat <- ldat$spliced
hist(log10(colSums(emat)),col='wheat',xlab='cell size')

# this dataset has already been pre-filtered, but this is where one woudl do some filtering

emat <- emat[,colSums(emat)>=1e3]

pagoda2处理 pagoda2用于生成细胞嵌入、细胞聚类以及更精确的细胞距离矩阵。您也可以使用其他工具生成这些工具，如 Seurat等。

创建pagoda2对象，调整方差：

library(pagoda2)
r <- Pagoda2$new(emat,modelType='plain',trim=10,log.scale=T)
#2600 cells, 7301 genes; normalizing ... using plain model winsorizing ... log scale ... done.
r$adjustVariance(plot=T,do.par=T,gam.k=10)
#calculating variance fit ... using gam 137 overdispersed genes ... 137 persisting ... done.

运行基本分析步骤以生成细胞嵌入和聚类，并可视化：

r$calculatePcaReduction(nPcs=100,n.odgenes=3e3,maxit=300)

running PCA using 3000 OD genes .. Loading required package: irlba .. done

r$makeKnnGraph(k=30,type='PCA',center=T,distance='cosine');

r$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=T)

绘制嵌入、集群标记（左）和"Xist"表达图（将男性和女性分开展示）

par(mfrow=c(1,2))
r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=10,shuffle.colors=F,mark.cluster.cex=1,alpha=0.3,main='cell clusters')
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$depth,main='depth')  

treating colors as a gradient with zlim: 1000.9 2939

速率估计

准备矩阵和聚类数据：

emat <- ldat$spliced; nmat <- ldat$unspliced; # disregarding spanning reads, as there are too few of them
emat <- emat[,rownames(r$counts)]; nmat <- nmat[,rownames(r$counts)]; # restrict to cells that passed p2 filter
# take cluster labels

cluster.label <- r$clusters$PCA$multilevel # take the cluster factor that was calculated by p2
cell.colors <- pagoda2:::fac2col(cluster.label)

# take embedding form p2

emb <- r$embeddings$PCA$tSNE

除了聚类和 t-SNE 嵌入，从 p2 处理，我们将采取细胞距离，优于默认的R 通常会使用全转录体相关距离。

cell.dist <- as.dist(1-armaCor(t(r$reductions$PCA)))

基于最低平均表达量（至少在其中一个cluster中）筛选基因，输出生成的有效基因总数：

emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.2)
nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
length(intersect(rownames(emat),rownames(nmat)))

[1] 2541

估计RNA速率：

fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=fit.quantile)

可视化 t-SNE 嵌入上的速率：

show.velocity.on.embedding.cor(emb,rvel.cd,n=200,scale='sqrt',cell.colors=ac(cell.colors,alpha=0.5),cex=0.8,arrow.scale=3,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1)

特定基因可视化：

gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 25,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,old.fit=rvel.cd,do.par=T)

调整视图

调整kCells，这会给我们一个更理想化的图像视图，差异肉眼可见：

gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 100,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,do.par=T)

文中链接

[1]

pagoda2: https://github.com/hms-dbmi/pagoda2

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注